Local Coulomb versus Global Failure Criterion for Granular Packings 
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Contacts at the Coulomb threshold are unstable to tangential perturbations and thus contribute 
to failure at the microscopic level. How is such a local property related to global failure, beyond the 
effective picture given by a Mohr-Coulomb type failure criterion? Here, we use a simulated bed of 
frictional disks slowly tilted under the action of gravity to investigate the link between the avalanche 
process and a global generalized isostaticity criterion. The avalanche starts when the packing as a 
whole is still stable according to this criterion, underlining the role of large heterogeneities in the 
destabilizing process: the clusters of particles with fully mobilized contacts concentrate local failure. 
We demonstrate that these clusters, at odds with the pile as a whole, are also globally marginal with 
respect to generalized isostaticity. More precisely, we observe how the condition of their stability 
from a local mechanical proprety progressively builds up to the generalized isostaticity criterion as 
they grow in size and eventually span the whole system when approaching the avalanche. 
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Understanding the failure of granular packings is of 
tremendous importance from both practical and theoret- 
ical aspects. Practically, avalanches are clearly of spe- 
cial interest for industrial and natural processes. From a 
more fundamental point of view, the mechanical rigidity 
of granular packings is related to the recently explored 
field of rheology close to dynamical arrest [H-IH , as well as 
to the nature of the jamming transition for frictional par- 
ticles 043 ■ Despite many studies both from a continuum 
and a microscopic point of view (see, for example, (8- 17|). 
the mechanisms of failure in frictional granular media are 
still unclear. 

Macroscopically, the application of the well known 
Coulomb criterion [l8| requires the knwowledge of an ef- 
fective friction coefficient, which remains out of reach of 
most recent developments. From a microscopic perspec- 
tive, Maxwell derived a global stability criterion based 
on counting the number of independent contact force 
components, which has to exceed the number of de- 
grees of freedom for a packing to be mechanically sta- 
ble [l9[ . Recently, this isostaticity criterion has been gen- 
eralized for frictional packings by including contacts ex- 
actly at the Coulomb threshold in the above counting ar- 
gument 2^, 21 1 . Such fully mobilized contacts are prone 
to tangential slipping, and it was indeed shown by Staron 



et al. [13, [15( that they play a key role in the destabiliza- 
tion process. A frictional packing of N particles with 
mean contact number z and mean number of fully mo- 
bilized contacts per particle n m has only Ndz/2 — Nn m 
independent force components, due to the additional re- 
strictions on the tangential forces. For this packing to 
be stable, this number has to exceed Nd(d + l)/2, the 
number of degrees of freedom and the the generalized 
isostaticity criterion in d dimensions reads: 
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where z^ so = d + 1 is the isostatic value for fj, = oo fric- 



FIG. 1: a - Generalized isostaticity phase diagram in 2D. In 
red is the line of marginal stability (Eq. [TJ. Frictional pack- 
ings under isotropic compression unjam - when compression 
is released (schematic green arrow) - at different positions 
close to this line (crosses) depending on friction \i [2(1 
Does a pile slowly inclined under gravity towards avalanche 
follow the hypothesized blue line? b - The system during 
the avalanche, at an inclination of 24° , for the details of the 
coloring scheme please refer to Fig. [3] 



tional packings [20| . It can be represented in a (z,n m ) 
phase diagram (see Figure [lj where a line of marginal 
stability divides the stable from the unstable regions of 
phase space. Recent molecular dynamics simulations us- 
ing an isotropic compression protocol (2p| have shown 
that frictional packings unjam (in the sense that p — > 0) 
close to the generalized isostaticity line. The final state is 
characterized by z(/jl) and n m (fi) and the linear response 
properties of these packings suggest that it is the distance 
to the line of marginal stability, 6z gen = z — _zf s e ™ , which 
controls stability [22j. 

These promising results, regarding both the global sta- 
bility criterion and the details of the microscopic mecha- 
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nism, raise two important questions. First, does the sta- 
bility criterion 5z 9en > remain valid in more realistic 
situations, which necessarily involve finite displacements 
and anisotropy? More specifically does a granular layer 
inclined under gravity follow the hypothetical trajectory 
plotted in blue in the (z, n m ) space on Figure[Tp Second, 
what is the link between the microscopic role of fully mo- 
bilized contacts in the failure of the system and the global 
generalized isostaticity criterion? 

In this paper, we explore these two issues using simula- 
tion data obtained by Deboeuf et al. Having precisely 
isolated the 'avalanches' from the 'quiet periods', our first 
immediate observation is that the pile destabilizes when 
it is still stable according to the generalized isostaticity 
criterion. This result prompted us to investigate the mi- 
croscopic role of the fully mobilized (or critical) contacts. 
These contacts form elongated clusters, the size of which 
obeys critical scaling, with a characteristic length which 
approaches system spanning size near the avalanche on- 
set. These observations, which are consistent with those 
by Staron et al. J13| , allow us to investigate the local and 
global stability properties of these clusters. Local failure 
as measured by the number of lost contacts correlates 
both spatially and temporally with being part of these 
clusters. Applying then the counting argument at the 
scale of these clusters, we find that their stability con- 
dition builds up progressively from a local mechanical 
property to the generalized isostaticity criterion as they 
grow in size. Hence the avalanche can be related to a sub- 
set of the packing becoming marginally stable according 
to the global criterion, while growing up to the system 
size. Note however that the above scenario is to be un- 
derstood as an averaged picture, since these clsuters are 
permcnantly destabilized and renewed. 

The system we study here is a simulation of 2d pack- 
ings of grains under gravity. The simulations where per- 
formed by Deboeuf et al. [H using the contact dynam- 
ics code developed by Staron (23j| . which assumes per- 
fectly rigid grains interacting at contacts through a hard 
core repulsion and a Coulomb friction law: the tangen- 
tial force at contact, ft, is related to the normal force /„ 
by the inequality \ft\ < where fi = 0.5 is the fric- 
tion coefficient. Beyond the fact that contact dynamics 
treats them as strictly nonsmooth, these contact laws 
do not differ from those more commonly used in dis- 
crete simulations [24| . The system consists of 4000 circu- 
lar grains with diameter uniformly distributed between 
[dmin, d-max] in a way to ensure 20% polydispersity. The 
length of the box is about 120d and the height is about 
35<i, where d is the mean diameter of the grains. Initially 
the grains free fall in the box set horizontally. The box 
is then tilted quasistatically to the desired inclination 9. 
Deboeuf et al. [f| have run several histories of inclination 
- including several back and forth oscillations - to inves- 
tigate the stress anisotropy. In the present case, we have 
selected the final part of the pile history, starting with a 
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FIG. 2: a - Distribution of the kinetic energy of the system 
averaged over all frames and all runs, b - Avalanche criterion: 
Probability density of the fraction of contacts lost between 
frames (log- binning). The criterion at 2% of failing contacts 
is shown in red. c - Distance from generalized isostaticity for 
a sample run; the run is divided into quiet and avalanching 
periods (in orange), d - Probability density of the position of 
the system in the generalized isostaticity diagram, overlaid by 
the trace of a sample run. The quiet period is shown in black, 
while the avalanching period is in red. The system clearly 
loses stability before crossing the marginal stability line, and 
the dynamics during the quiet period is dominated by the 
fully mobilized contacts. 



horizontal pile with 9q = and tilting it in the direction 
of positive 9 until 9 w 30 . We use 50 independent runs, 
with different initial conditions and our results are usu- 
ally shown as an average over these runs. For each run, 
particle positions, contacts and forces were stored in suc- 
cessive frames separated by 50 computational time steps, 
which corresponds to an angle variation of 59 - 0.05 °. 

Let us first turn to the identification of the avalanche 
periods. We are interested in the properties of the pack- 
ing just before failure occurs, i.e. for 9 — > 9 m , where 
9 m , the failure angle is, known to show considerable sta- 
tistical variance [5|. Obtaining a good criterion to pre- 
cisely identify the start of an avalanche is not straight- 
forward: one possibility with limited resolution is to con- 
sider the probability distribution of the kinetic energy 
(cither translational and rotational) and identify as a 
threshold the value above which the distribution escapes 
the power law reported in the quiet region (see the arrow 
in Figure^). A more explicit signature of failure is given 
by considering the probability density of the fraction of 
contacts lost between two subsequent snapshots. As ob- 
served in Figure [2jb) , the data clearly fall into two dis- 
tinct sub-populations (note the logarithmic scale on the 
x-axis), which we associate with the quiet period (left) 
and the avalanching period (right). The good separation 
between the two populations allows us to apply a cut- 
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off of a maximum of 2% of contacts lost between frames 
to determine the end of the quiet period (in red). In 
Figure [5Jc) , this criterion has been applied to a sample 
configuration, and it is clearly able to capture the loca- 
tion of significant events in the packing. We have checked 
the robustness of this criterion by varying the cutoff be- 
tween 1% and 4%, which typically changes the location 
of the onset of the avalanche by 1-2 frames, correspond- 
ing to an error in 9 m of 0.05° — 0.1°. We have observed 
that the position of the avalanche onset is very variable, 
and that there may be several avalanches of various sizes 
within a single run. In the following all ensemble aver- 
ages are performed as a function of 9 m — 9, where 9 m is 
the location of the next avalanche onset. 

Figure [5Ji shows the probability density of the system 
in the parameter space (z, n m ), together with the trace of 
one sample run. The black (respectively the red) symbols 
correspond to the quiet (resp. avalanche) period. Two 
features are apparent. First, the quiet period and the 
avalanche occupy very well separated regions of the pa- 
rameter space and the transition between them is abrupt. 
During the quiet period, the mean contact number z is 
quenched, and all of the dynamics is due to fully mobi- 
lized contacts appearing and disappearing in the packing 
(corresponding to an up-and-down motion in the graph). 
This is confirmed by our observation that throughout the 
pre-avalanche period, particles move typically less than 
1% of their diameter and remain "caged". When the 
system enters the avalanche period, the structure of the 
packing breaks down and the trace moves along a diag- 
onal. Second, the system crosses the marginal stability 
line only after the start of the avalanche. This means 
that the generalized isostaticity criterion, (Eq. [I}, is a 
necessary criterion for stability, but not a sufficient one. 

Generalized isostaticity is a global, mean-field type cri- 
terion, and it is likely to overestimate the stability of the 
packing if the system behaves in an heterogeneous way. 
We have checked that such heterogeneities are not related 
to the geometry. If we exclude the top and bottom sin- 
gle layers of particles, which cannot easily be integrated 
into the counting, no part of the system is critical at 
the start of the avalanche (note that we only used the 
bulk of the system for Figure [2J1). We also have cut the 
system both into lanes, and several other patterns and 
found that the same general conclusion holds: within a 
scale of a few frames, the avalanche onset measured by 
this method is the same for all sub-parts of the system 
ensuring that the avalanche is truly a global property of 
the system (see also Figure HJd) . Note that for a granu- 
lar bed much deeper than our relatively shallow layer of 
particles, one may observe a stronger dependance of n m 
on depth than in the present case and that the system 
could then eventually separate. 

Inspired by Staron's work (l3j |. we identify the crit- 
ical contacts as a natural candidate for the source of 
the heterogeneity suggested by the above observations. 
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FIG. 3: (Color online) Top: Contact properties close to the 
avalanche onset, 8 m — ~ 0.3°; clusters of particles with at 
least one critical contact (in green) are in black. Particles with 
failing contacts belonging to a cluster are red, failure of the 
other particles is in blue. The red arrow indcates the direction 
of gravity - note the directionality of the clusters. Bottom: 
Spatial (left) and temporal (right) correlation between clus- 
ters and failure during the quiet period. Left: Fraction of 
failure which occurs in a cluster, N[ /N c as a function of frac- 
tion of the system in a cluster, N c /N. The dashed line is 
the limit where failure and clusters are uncorrelated. Right: 
Temporal correlation between the maximum cluster size at 
time t and the number of particles that fail between t and 
t + t; the units are frames each corresponding to A9 = 0.1°. 



Since the contact number is quenched during the quiet 
period and all changes in Sz gen are borne by the criti- 
cal contacts, while the remainder of the system does not 
contribute, we can divide the system into two subpopu- 
lations: those particles with at least one fully mobilized 
contact, and those without. Figure(3J top, illustrates the 
critical contacts close to the onset of the avalanche. It be- 
comes apparent that the particles with critical contacts 
organize in rapidly fluctuating clusters. These clusters 
are rather anisotropic and their preferential direction of 
alignment makes an angle of roughly 60° with respect to 
gravity. This is in agreement to the Mohr- Coulomb fail- 
ure criterion and its extension in the context of elasto- 
plasticity [25|, which would predict a failure angle of 
1/2 (ir/2 — atan(/Lt)) ~ 30°, that is along the slip plane 
at the critical contacts, normal to the contact orienta- 
tion. 

We now turn to the relation between these clusters and 
local failure, as measured by the contacts lost between 
two consecutive frames. We find both spatial and tempo- 



FIG. 4: Probability density of the position of the clusters in 
the generalized isosaticity diagram, before the avalanche in 
black, and during the avalanche in red. The clusters straddle 
the marginal stability line during the quiet period. The z- 
distribution of the remainder of the system is shown below - 
note that since n m = for this subset, we have Zi ao = 3. 

ral correlations between clusters and failing contacts. As 
illustrated in Figure [31 particles within clusters are more 
likely to have failing contacts. The particles with failing 
contacts which belong (don't belong) to a cluster are il- 
lustrated in red (blue) in the top figure, and even though 
less than half the system is part of a cluster, most of the 
failure occurs in clusters. The bottom-left figure confirms 
that the fraction of particles which fail and belonged to a 
cluster in the previous frame as a function of the fraction 
of the system which is part of a cluster, is much higher 
than if both events were uncorrelated (dashed line). Sim- 
ilarly, we track the size S m of the largest cluster in the 
system and correlate it with the number of contacts N* 
that fail in the following frames using 

r( , = J2(S m (k)-(S m ))(Nf(k+r)-(Nf)) 
[T> {Var{S m )Var(Nf)y/ 2 ' 1 ' 

where the limits on the sum are chosen such that both 
k and fc + r fall within the quiet period. The mean and 
the variance are computed over the same set of frames. 
As shown in Figure [3] bottom-right, both quantities are 
correlated during the whole quiet period, but the maxi- 
mum correlation happens at r = 0: just after the cluster 
size peaks, the failure rate is enhanced by about 40%. 
The slow decay for r > is due to the persistence of 
clusters, since most of critical contacts tend to belong to 
particle pairs which repeatedly form and loose contact. 
In contrast, during the avalanche, both the correlations 
between clusters and failure and the directionality of the 
clusters are lost (see Figure [J) - a signature of the break- 
down of the structure of the system. 

The clusters of particles with at least one critical con- 
tact hence appear as natural candidates for being the 
seeds of the destabilization process. However, for sta- 
bility reasons, this is far from trivial. Let us recall here 



FIG. 5: (Color online): Cluster size distribution p(s) at dif- 
ferent distances from the avalanche onset, ensemble averaged 
over an angular interval of AO — 0.5°. The arrow is in the 
direction of approaching the avalanche. From left to right, 
(9 — 6 m ) - values: —14.4° (green plusses), —4.6° (blue trian- 
gles), —0.75° (red squares), —0.25° (black dots); and during 
the avalanche (purple diamonds). Lines are fits to Eq. O 
and the dashed line indicates a slope of —2.5. Insets: (a) £, 
(b) maximum cluster size S ma x as a function of 8 — 9 m , the 
distance to the avalanche onset. 



that for one particle to be stable, it must satisfy the local 
stability criterion z — 2n m /d >l/2(d+l), which states 
that the dz — 2n m force components need to be able to 
constrain the l/2d(d + 1) degrees of freedom of the par- 
ticle. The particles which belong to these clusters have 
by construction a minimal n m = 0.5 (at least one critical 
contact shared by two particles) . The average number of 
contacts for the whole pile is (z) ~ 3.5. Thus these par- 
ticles easily satisfy the local stability criterion. However 
when such individually stable particles aggregate into 
clusters, one expects that above some mesoscopic size, 
the cluster will need to satisfy the global criterion ([1}. 
Assuming for the moment that the average contact num- 
ber within the clusters is also close to (z) ~ 3.5, one ob- 
tains that these clusters are prone to be marginally stable 
as confirmed in Figure 01 where one can see that the clus- 
ters straddle the marginal stability line during the quiet 
period (in black) and hence arc highly unstable to per- 
turbations. And indeed one observes visually that they 
are very intermittent, permenantly loosing and gaining 
particles, vanishing and reforming in contrast with the 
remainder of the system is always far removed from iso- 
staticity during the quiet period. Hence, it is not clear 
that the clusters can grow up to a system spanning size. 

As a matter of fact, this complex dynamics leads to 
an interesting critical feature for the the cluster size dis- 
tribution. As shown in Figure [5] it develops larger and 
larger tails when approaching the avalanche onset. These 
distributions are indeed well fitted by a power law with 
an exponential cut off, the characteristic lengthscale of 
which increases sharply when approaching the avalanche 
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FIG. 6: Cluster size effects on stability: a - Distribution of Sz 
for different cluster sizes, with binning that reflects the fun- 
damental discreteness of the distributions. The curves have 
been shifted for clarity, b - In black (squares) is {Sz) for differ- 
ent cluster sizes and the global stability criterion (dash-dot). 
The blue dots are the minimal Sz obtained for each cluster 
size; the dashed line is the marginal stability criterion for a 
single grain. 



onset: 



P(s) 



(3) 



where s is the size of the cluster in particles. We 
have estimated 5 = 2.5 ± 0.25 from the distribution 
just above onset (see dashed line) and the inset (a) dis- 
plays the cluster scale £ as a function of the distance 
to avalanche onset. We observe a sharp upturn close 
to the avalanche, especially in the last 10 frames. Sim- 
ilarly, the ensemble-averaged size of the largest cluster 
also grows and presents an upturn just before the onset 
of the avalanche (see inset (b) of Figure [5]) . Hence, de- 
spite their intrinsic marginal stability, and provided that 
the clusters are indeed quasi unidirectional, which was 
also established in (l5| . the maximal size recorded here 
corresponds to clusters reaching system-spanning size at 
the onset of the avalanche. 

In the above discussion, we have assumed that the av- 
erage contact number inside the clusters is roughly equal 
to the one for the whole pile. However given the in- 
crease of the exponential cut-off when approaching the 
avalanche, one suspects that the average contact num- 
ber and to a lesser extent the average number of critical 



contacts actually depend on the cluster size. To elu- 
cidate this last point, we compute the distributions of 
Sz = z — 2n m /d within clusters of a given size N (see 
figure [(Ji). Apart from discretization effects at small size 
and the shrinking of the width with ~ TV^ 1 / 2 , the distri- 
butions progressively shift from negative Sz to a slightly 
positive one. Note also that the minimal value of Sz 
goes from the marginal local criterion for the smallest 
one (all studied clusters must be composed of locally sta- 
ble particles) to close to the global stability criterion for 
the largest one (figure 0b). Altogether these clusters ap- 
pear to be dynamically selected according to their stabil- 
ity, which from a local mechanical proprety progressively 
builds up self consistently with their size towards the 
generalized isostaticity criterion. 

The picture that emerges from this analysis is the fol- 
lowing. The idea that the whole system is isostatic at 
the avalanche onset is too simple because it ignores the 
anisotropy and inhomogeneity of the pile, which trans- 
poses to the critical contacts. Such important effects 
of the anisotropy have been strongly emphasized be- 
fore [1, HH and it is clearly a too strong assumption 
to elude them [2f|. This suggests an even deeper his- 
tory dependence of frictional piles, which calls for a re- 
fined description of the texture beyond the introduction 
of the number of critical contacts. For the isotropically 
compressed packings of [13, HH, it was also noted that 
only extremely slowly equilibrated packings unjam at 
S?9 en = 0; for those packings, the critical contacts are 

The appearance of 



27 



indeed randomly distributed 
system-spanning marginally stable clusters is an intrigu- 
ing mechanism of unjamming; one which cannot exist in 
frictionlcss systems if the packing structure and hence the 
local contact numbers are to remain homogeneous. This 
suggests that the frictional and the frictionless jamming 
transition may be more different than expected. Fur- 
ther insight in that matter could be gained from new ex- 
periments with two dimensional packings of photoelastic 
discs; in particular, the generalized isostaticity diagram 
could have direct relevance to the cyclic shear experi- 
ments of (H. 
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